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It is well known that the Black-Scholes-Merton model suffers from several deficiencies. Jump- 
diffusion and Levy models have been widely used to partially alleviate some of the biases inherent 
in this classical model. Unfortunately, the resulting pricing problem requires solving a more 
difficult partial-integro differential equation (PIDE) and although several approaches for solving 
the PIDE have been suggested in the literature, none are entirely satisfactory. All treat the 
integral and diffusive terms asymmetrically and are difficult to extend to higher dimensions. We 
present a new, efficient algorithm, based on transform methods, which symmetrically treats the 
diffusive and integrals terms, is applicable to a wide class of path-dependent options (such as 
Bermudan, barrier, and shout options) and options on multiple assets, and naturally extends 

>• to regime-switching Levy models. We present a concise study of the precision and convergence 
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properties of our algorithm for several classes of options and Levy models and demonstrate that 



the algorithm is second-order in space and first-order in time for path-dependent options. 
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1. Introduction 



The seminal works of Black and Scholes (1973) and Merton (1973) (BSM model) revolution- 



ized our understanding of financial contracts with embedded options. Based on the assumption 
that stock prices follow a geometric Brownian motion, i.e. stock returns have a log-normal dis- 
tribution, they demonstrated that a replicating strategy reduced the pricing problem to solving 
a partial differential equation (PDE) which is independent of the return of the asset. Today, 
option traders quote prices in terms of implied volatilities induced by matching market prices 
with those of the BSM model; however, it is well known that the BSM model suffers from sev- 
eral deficiencies rendering it inconsistent with market price behaviour. These inconsistencies 
manifest themselves in, for example, the observed implied volatility smile (or skew) and term 
structure. Various lines of research aim to remove these pricing biases by focusing on disparate 
extensions. One line of research seeks to introduce state dependence resulting in correlations 



between prices and volatility levels (see e.g. Derman and Kani (1994) and Duprie (1994) for 
nonparametric models and Cox and Ross (1976)| Ingersoll (1997) and |Rady (1997) for para- 



metric models). Another line of research elevates volatility to a stochastic variable (e.g. Heston 



(1993) I, or assumes volatility undergoes regime changes (e.g. Naik (1993) ). A third line of 



research focuses on introducing jumps into the prices process itself (e.g. iMerton (1976) and 



Madan and Seneta (1990) ) while maintaining time homogeneity. All of these directions are able 



to correct for different aspects of the implied volatility surface and have their own unique set of 
advantages and disadvantages. 

In this paper, we focus on pricing options where the underlying index, or indices, are driven by 
Levy processes both with and without regime changes. This combines two of the three modeling 
directions and we succeed in developing an efficient method of pricing for a wide class of options. 
In all, there are four main purposes for this paper: first, to introduce a new numerical method 
based on the Fourier transform of the pricing partial integro-differential equation (PIDE); second, 
to study the algorithmic performance for various European and path-dependent options; third, to 
extend the method for multi-asset path-dependent contingent claims; and, fourth, to incorporate 
regime changes. 

The option pricing problem under the BSM model can be reduced to solving a second-order 
parabolic PDE with the independent variables being time and stock price. By changing terminal 
or boundary conditions, or imposing early exercise constraints, the PDE can be used to price 
a variety of options. Under jump models, a PIDE with a non-local integral term must now 
be solved. A quick review of exponential Levy models and the pricing PIDE is provided for 
completeness in section [2j An assortment of finite difference methods for solving these PIDEs 



have been proposed in literature, see e.g. Andersen and Andreasen (2000) Briani, Natalini, and 



Russo (2004)1 |Cont and Tankov (2004)] and |d'Halluin, Forsyth, and Vetzal (2005) [ Although the 
methods are quite diverse, they all treat the integral and diffusion terms of the PIDE separately. 
Invariably, the integral term is evaluated explicitly in order to avoid solving a dense system of 
linear equations. In addition, the Fast Fourier Transform (FFT) algorithm is employed to speed 
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up the computation of the integral term (which can be regarded as a convolution) and/or its 
inverse. Unfortunately, these methods require several approximations such as: 

• in infinite activity processes, small jumps are approximated by a diffusion and incorporated 
into the diffusion term; 

• the integral term must be localized to the bounded domain of the diffusion term, i.e. large 
jumps are truncated; 

• the option price behaviour outside the solution domain must be assumed; 

• the separate treatment of diffusion and integral components requires that function values 
are interpolated and extrapolated between the diffusion and integral grids in order to 
compute the convolution term. 

These factors together make finite difference methods for option pricing under jump models 
quite complex, and potentially prone to accuracy and stability problems, especially for path 
dependent claims. 

In section [3j we present a new Fourier Space Time-stepping (FST) algorithm. This method 
avoids the problems associated with finite difference methods by transforming the PIDE into 
Fourier space. One of the advantages of working directly in Fourier space is that the character- 
istic exponent of a independent increment stochastic process can be factored out of the Fourier 
transform of the PIDE. Consequently, the Fourier transform can be applied to the PIDE to 
obtain a linear system of easily solvable ordinary differential equations (ODE). Furthermore, 
the characteristic exponent is available, through the Levy-Khintchine formula, in closed form 
for all independent increment processes. This makes the FST method quite flexible and generic 
- contingent claims on any exponential-Levy stock price processes can then be priced with no 
additional modifications to the algorithm. The FST naturally leads to a symmetric treatment 
of the diffusion and jump terms and avoids any explicit assumptions on the option price outside 
of a truncated domain. 

Since the FST method provides exact pricing results between monitoring times, it is signif- 
icantly more efficient and accurate when compared with finite-difference methods for valuing 
Bermudan options. Furthermore, the method allows prices from one monitoring time to be 
projected back to a second monitoring time in one step of the algorithm. Contrastingly, finite- 
difference schemes will require time-stepping between monitoring dates resulting in further pric- 
ing biases and speed reduction. 

For path independent options, prices for a range of spots can be obtained in a single time 
step. For exotic, path dependent options, we demonstrate how the FST method can handle 
barrier and Bermudan (American) styled clauses. The closed form expression for the Fourier 
transform of the option payoff is not required, making the FST method easily applicable to 
options with non-standard payoffs. Since the FST method requires two FFTs per time step, its 
computational complexity is 0{K N logiV), where N is the number of spatial gridpoints and 
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K is the number of time steps. Through numerical experiments, we establish that the order of 
convergence of the method is two in space and one in time for path dependent options. 

In section |4j we generalize the method to the multi-asset case and apply it to the specific 
example of spread options and comment on its use in catastrophe option pricing. In typical 
markets, jump models alone cannot match the implied volatility skew for longer maturities; 
however, the observed market behaviour can be captured by incorporating regime switches. This 
motivates us to include one more generalization and we introduce a non-stationary extension 
of the multi-dimensional Levy processes using regime changes. The regime changes are induced 
through a homogenous continuous time Markov chain. This allows the index(es) to exhibit 
stochastic volatility and/or stochastic correlation behaviour which can be important in longer 
term options. Stochastic correlation has received little attention in the literature to date; however 
our modeling and pricing framework easily handles this feature. 

We conclude this paper by discussing the possible applications, improvements and extensions 
of the FST algorithm. 



2. Option Pricing with Differential Equations 

In this section, we review the differential equation approach to option pricing with exponential 
Levy processes. For a modern treatment of this subject the interested reader is referred to the 



monograph by Cont and Tankov (2004) and to Sato (1999) for further mathematical background. 



Let V(t, S(t)) denote the price at time t of an option, written on an underlying price index S(t), 
with a T-maturity payoff of (p(S(T)). It is well known that, in an arbitrage-free and frictionless 
market, the value of the option is the discounted expectation under a, not necessarily unique, 



risk-neutral measure Q (see Harrison and Pliska (1981) ). Explicitly. 



V(t,S(t))=E?[e-^ T - t ^(S(T))\ , (1) 

where the expectation is taken with respect to the information, or filtration, Tt available at 
time t. Here and in the remainder of this article, we assume that the risk-free interest rate r 
is constant. When the underlying index follows a diffusion process, the risk-neutral measure 
is indeed unique; however, in the more interesting case of exponential Levy models, many 
equivalent risk-neutral measures exist. Nonetheless, we take the view that a trader is using such 
a model to price derivative instruments and therefore is modeling directly under a particular 
risk-neutral measure - possibly induced through a calibration procedure. 

A dual and equivalent specification of the value function is its associated PIDE formulation. 
These two specifications are connected by noting that the discount-adjusted and log-transformed 
price process v(t, X(t)) := e r V(t, S(0) e x ^) is a martingale under the measure Q. Conse- 
quently, the associated drift term of its defining SDE is identically zero. If the underlying index 
follows an exponential Levy process, then the price process can be written as S(t) = S{0) e x ^ 
where X(t) is a Levy process with characteristic triplet (7, cr, z/). In this case, the process X(t) 
admits the following canonical Levy-Ito decomposition into its diffusion and jump components 
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(see |Sato (1999)] ): 

X(t) = ~ft + aW{t) + J l {t) + lim J e {t) 



e\0 



J\t) 



y\i{dy x ds) , 



o J \v i 

t 



Je<|y|<l 



y [//(<iy x ds) — v(dy x (is) 



(2) 
(3) 
(4) 



Here Wt is a standard Brownian motion, fJ,(dy x ds) is a Poisson random measure counting the 
number of jumps of size y occurring at time s, and v{dy x ds) = v{dy) ds is its compensator. 
Note that J (t) (J € (t)) carries the interpretation of large (small) jumps. If the model has finite 
activity (/]Ry{o}(M A 1) v[dy) < +oo), then there is no need to truncate small jumps and they 
can be lumped together with large jumps. We, however, choose to leave the decomposition 
general. Given this modeling assumption, v must satisfy the PIDE 



(d t + C)v 
v(T,x) 



0, 



(5) 



where C is the infinitesimal generator of the Levy process and acts on twice-differentiable func- 
tions f(x) as follows 



Cf(x) 



lim 

t\o 



E[f(x + X(t))]-f(x) 



jdxf + \a 2 d xx f + / [f(x + y)- f{x) - y 1 {M<1} d x f(x)] u{dy). (6) 
iR/{0} 



By enforcing the risk-neutrality condition, the drift is uniquely determined once the volatility 
and Levy density are specified. In particular, 7 satisfies 



En 



1 



(7) 



where ^Sf(uj) denotes the characteristic exponent of the Levy process and is provided explicitly 



by the Levy-Khintchine formula (see Sato (1999)) 



:=lnE^[e 



!U)X(l)l 



?7W 



a lUJ V 



1 



R/{0} 



^y^{\ y \<i})v{dy). 



(8) 



Within this framework, the classical purely diffusive (BSM) model is recovered by setting the 
Levy density to zero. Furthermore, jump-diffusion models, in which the log-stock price contains 
a diffusive component together with jumps occurring at Poisson times, are recovered by setting 
v(dy) = A fy(y) dy where A is the activity rate of the Poisson process and fy(y) is the probability 
density of the jumps. In this case, the process X(t) can be written in terms of a Q-standard 
Brownian motion W(t), a Poisson process N(t) with activity rate A, and i.i.d. random variables 
Yi, representing the jumps at Poisson times U, as follows: X(t) = 7 1 + a W(t) + J2n=i ^n- Two 
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Model 



Levy density v{dy) 



Characteristic Exponent ^(uj) 



BSM 

Merton JD 
Kou JD 

VG 

NIG 
CGMY 
Table 1 
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The Levy densities and characteristic exponent for various models. Here a = (3 = — — , 

-y/ (J, 2 +<J 2 /k „ _ -y/ /i 2 +(7 2 /K 



7 



and Kp(z) is the modified Bessel function of the second kind. 



widely used jump- diffusion models are the log-normal jump model due to Merton (1976) and 



the double exponential model due to Kou (2002) 



The Kou (2002) model assumes that the X(t) jumps are double exponentially distributed, 
with positive jumps (with probability p) of mean size r/ + and negative jumps (with probability 
1 — p) of mean size rj- . 

The characteristic exponents and Levy densities for these models are provided in Table [TJ 
These jump-diffusion models are popular not only because they perform well when calibrating 
to option prices, but also because they admit two semi-explicit closed form solutions. One 
form involves an infinite summation of BSM like prices (which can be safely truncated to a 
small number of terms), while the other form involves an inverse Fourier transformation. The 
interested reader is referred to the respective papers for details. 

More recently, pure jump models have become very popular across a number of markets in- 
cluding equity, interest rate and commodity markets. These models have been found to better 



fit implied volatility smiles than jump-diffusion models and are widely used in industry. Huang 



and Wu (2004j| carry out numerous statistical tests which demonstrate that models with in- 



finitesimal jumps outperform jump-diffusion models. Within this class, the jumps themselves 
occur infinitely often with most jumps being of infinitesimal size. Several breeds of pure jump 
models have been suggested in the literature and each has its own merits and drawbacks. Three 
very popular models are the Variance-Gamma (VG) model of Madan and Seneta (1990) and 



Madan, Carr, and Chang (1998) the CGMY extension of the VG model developed by Carr, 



Geman, Madan, and Yor (2002) and the normal inverse Gaussian (NIG) model popularized by 



Barndorff-Nielson (1997). The various Levy densities and characteristic functions are provided 
in Table [TJ 
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In the absence of jump components, the resulting PDE can be discretized using standard 
divided differences to approximate the first and second order derivatives. Whether the approxi- 
mation scheme is carried out explicitly or implicitly or through a weighted scheme, the resulting 
system is tri-diagonal and leads to very efficient numerical approximations. Unfortunately, when 
jumps are present, the integral term in ([5j must be approximated resulting in a dense matrix 
structure. Several methods for dealing with this issue have been presented in the literature with 
most relying on the explicit evaluation of approximations to the integral term in conjunction 
with an iterative refinement and possibly FFT speedup. 



Andersen and Andreasen (2000) propose an FFT-Alternating Direction Implicit (FFT-ADI) 
method which treats the diffusion and integral terms symmetrically over a full time step by 
splitting the time step into two half-steps; an explicit scheme is used on the first half-step 
and implicit scheme on the second half-step. The inversion of the dense matrix is performed 
efficiently by regarding the term as a convolution and utilizing the FFT algorithm. The fixed 



point iteration scheme of d'Halluin, Forsyth, and Vetzal (2005) treats the integral term explicitly 



while iterating to attain required error tolerance per time-step. In addition, Implicit-Explicit 
(IMEX) Runge-Kutta schemes have been applied in Briani, Natalini, and Russo (2004)] to solve 



the PIDEs. Although the FFT algorithm is frequently used to speed up the computation of the 
integral term, such schemes require careful mapping of function values between the diffusion 
and integral grids. We circumvent the problems posed by working in real space and instead opt 
to solve the problem directly in Fourier space as explained in the next section. 

3. Fourier Space Time-stepping 

Fourier and Laplace transforms have been used extensively to solve PDEs, either by trans- 



forming the equation into an ODE or expressing the solution as an infinite series (see Strauss 



(1992) and Taylor (1997)). The aim of this section is to develop a Fourier transform method- 
ology for solving PIDEs of the form ([5]). The main advantage of transform methods is that the 
PIDE can be handled efficiently without the additional complexities associated with the integral 
term. Additionally, the algorithm is applicable to any independent increment stock price model 
which admits a closed form characteristic function. Furthermore, we extend this approach to the 
valuation of path dependent options, such as barrier, American, and shout options, and discuss 
the convergence of these numerical schemes. 

3.1. Transforming the PIDE 

A Pseudo Differential Operator (PDO) extends the notion of a differential operator and is 
widely used to solve differential equations. The essential idea is that a differential operator with 
constant coefficients can be represented as a composition of a Fourier transform, multiplication 
by a polynomial function, and an inverse Fourier transform. Only a few fundamental facts from 
PDO theory are required to derive our numerical method. The interested reader is referred 
to Boyarchenko and Levendorskii (2002) who discuss the PDO theory in the context of option 



pricing. For a more thorough treatment of the subject see Taylor (1997) 



8 



K. Jackson, S. Jaimungal and V. Surkov 



A function in the space domain fix) can be transformed to a function in the frequency domain 
J~\f\(oj) (where uj is given in radians per second) and vice- versa using the Fourier transform: 

/oo -i roo 

f{x)e-^dx and F~ x \f\{x) x= — I f[u:)e^dw. 
-oo •^ 7r J — oo 

The Fourier transform is a linear operator that maps spatial derivatives d x into multiplications 
in the frequency domain: 

F[i%f\ (t,w)=iwf[dr 1 f] (t,w) = ■■■ = iiwmm.w). 

Consequently, applying the Fourier transform to the infinitesimal generator C of X{t), defined 
by equation allows the characteristic exponent of X(t) to be factored out: 

F[Cv][t,u) = |*7"-^+/ [e^ x -l-iuyt M<1} )u{dy)\r[v)(t,uj) 
{ 2 Jb./{o} J 

= *(u)F[v](t,u). (9) 
Furthermore, taking the Fourier transform of both sides of the PIDE §5§ leads to 

f dtr[v](t,u) + *(u)r[v](t,u) =o, 



The PIDE is therefore transformed into a one-parameter family of ODEs (10) parameterized by 



u. Giving the value at time t 2 < T, the system is easily solved to find the value at time t\ <t 2 : 

F[v](t 1 ,u;) = F[v}(t 2 ,u J )-e&- t ^\ (11) 
Taking the inverse transform leads to the final result 

v(tx,x) = T- x \T\v\{t 2 ,u) • e fe-*0*(«)} (x). (12) 
3.2. Direct Transform Method 



Alternatively, it is possible to derive (12) directly from the expectation representation of 



prices. Recall that v is a Q martingale; consequently, 

v(h,X(ti)) = E®[v(t 2 ,X(t 2 ))] 

v(t 2 ,X(h) +x) fx{t 2 )-x{t x ){x) dx 

) 

v(t 2 ,X(h) +x) f X {t 2 -t x ){x) dx ■ 



-oo 
oo 



Here, fx(t)( x ) denotes the p.d.f. of the process X{t) and the third line follows from the inde- 
pendent increment property of the process X(t). Furthermore, ,F[/x(i)]( w ) = e t *(-" ; ) and since 
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a convolution in real space corresponds to multiplication in Fourier space, we have 

nv](ti,u) = F[v](t 2 ,uj) eto-WW . (13) 
3.3. FST Method 



Armed with the pseudo-differential operator solution (12 ), the numerical algorithm is straight- 



forward. For path-independent options the price is obtained in one step by directly applying 



equation (12) similar in spirit to Carr and Madan (1999) For path dependent options a time- 



stepping algorithm is used to apply boundary conditions, impose constraints, or optimize over 
a policy domain. 

Consider a partition of the time and truncated stock price domain 0, = [0,T] x [x m i n ,x max ) 
into a finite mesh of points {t n \n = 0, . . . , iV} x {x m \m = 0, . . . , M — 1}, where t n = nAt, x m = 
Xmin+mAx and At = T/N, Ax = (x max -x min )/M . The mesh boundary {x min , x max } is chosen 
large enough to capture the overall behaviour of the option value function, yet small enough 
to maintain the accuracy of the computed option price in the range of interest. Numerical 
experiments suggest that x m i n S [—4, —2] and x max £ [2, 4] works well for diffusion models, while 
Xmin £ [ — 5, — 3] and x max £ [3,5] is preferable for models with a dominant jump component. 
Note that x = log(5/5o). Analogously, consider a partitioning of the time and frequency domain 
£1 = [0, T] x [0, LU max ] into a finite mesh of points {t n \n = 0, . . . , N} x {uj m \m = 0, . . . , M/2}, 
where u> m = mAuj and Aw = 2u> max /M. We choose ui m ax = §Sa;> w ^ich is the Nyquist critical 
frequency. Note that v(t,x) is a real- valued function and thus F\v\ (t, — u) = T[v\ (t, u) . The 
Fourier transform for negative frequencies is not computed and therefore the frequency grid has 
half as many points as the spatial grid. 

Let v m := v{t n ,Xm) represent v(t,x) at the node points of the partition of $7, and let v m := 
ni w m) represent J~[v] (t, to) at the node points of the partition $7. The frequency domain prices 
are obtained from the spatial domain prices as follows: 

Af-i 



k=0 

M-l 



OL m £ 4 v k e 

k=0 



= a m FFJ[v n ](m). (14) 

Here, a m = e~ luJmXmin Ax and FFT[f n ](m) denotes the m-th component of the discrete Fourier 
transform (DFT) of the vector v n , which can be computed efficiently using the FFT algorithm. 
Similarly, the spatial domain prices can be computed from frequency domain prices via a discrete 
inverse transform 

u^ = FFT- 1 [a- 1 .v n ](m). (15) 



Combining these connections between frequency and spatial domains with the transformed 
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Real Space Fourier Space Real Space Fourier Space 

FFT 



yn 



yn 



Apply B.C. 



F[V n ] 



U[V T ' 



FFT 



■F[U[V n ]] 



Time-step 



Apply B.C., FFT 



T\y n ■ h b ] = F\y n \ <g> f[h b ] 



Time-step 



FFT" 1 

V n ~ l ^[V"- 1 ] 



yn-l 



FFT 



Figure 1. A schematic representation of the FST algorithm. In the left panel, the boundary 
conditions (such as optimal exercise or barrier breach) are applied in real space while the time 
step is performed in Fourier space. In the right panel, the refined algorithm for Barrier breach 
permits application of boundary conditions in Fourier space for particular cases. 



PIDE (12), a step backwards in time is computed by 



, ; n-l = FFT — 1 [qe — 1 - w**- 1 ] 

= FFT -1 [a -1 • v n ■ e* At ] 

= FFT -1 [aT 1 • a • FFTfu' 1 ] • e* At ] 

= FFT _1 [FFT[f n ] • e* A ']. 



(16) 



Notice that the coefficient a, which embeds information about the spatial boundary, cancels 
in the above equation and can be omitted during the numerical computation. Certain payoffs 
contain singularities in their Fourier transforms along the real axis. A simple shifting of to — > 
u) + ie avoids this problem, resulting in a slight modification of the time-stepping algorithm: 



„n-l 



FFT~ 1 [FFT[# n ] • e* At ] where v% = e eXk v% and *(w) = *(w + ie). 



3.3.1. European Options 

European options can be valued in a single time step, since (16) is a valid approximation for 
any At. In this case, given a payoff function (p(S), set N = 1, = (p(S(0)e Xm ), numerically 



invert to obtain via (14), and finally apply (16). This approach is similar to Carr 



and Madan (1999) , however, the explicit expression of the Fourier transformed option payoff is 
not required - clearly a great advantage for non-standard payoffs. Moreover, our approach is 
computationally more efficient when compared to spatial PIDE solution-based methods since it 
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does not require stepping in time. 

3.3.2. American Options 

American options can be priced using a finite difference method either by solving a linear 



complementarity problem (see Dewynne, Howison, and Wilmott (1993), Dempster and Hutton 



(1997) , Huang and Pang (1998) and Forsyth and Vetzal (2002) ) , or solving a free boundary value 



problem (see McKean (1965) Kim (1990) and Carr, Jarrow, and Myneni (1992)). Although 



the free boundary formulation for American options is an active area of research and can be 
potentially combined with the FST method, the linear complementarity formulation is easier to 
implement in the context of the FST method. Since the value of an American option is always 
greater than or equal to the terminal payoff, the idea is to continuously enforce the condition 
V(t, S) > V(T, S). Numerically, this is enforced when boundary conditions are applied, resulting 
in the following algorithm: 



v 



n-1 



max {FFT _1 [FFT[f n ] ■ e q ' At ],v N } . (17) 



where the max(-, •) is taken componentwise. There is no convenient representation of the max(-, •) 
operator in Fourier space; consequently, it is necessary to switch between real and Fourier spaces 
at each time-step. Schematically, the algorithm is presented in the left panel of Figure [T] 

3.3.3. Barrier Options 

The numerical algorithm for barrier options is similar to that of American options and also in- 
volves enforcement of constraints. Here, we discuss the up-and-out barrier option case; however, 
the results can be extended to other barrier option styles. In spatial coordinates, the barrier 
boundary condition forces 

V(t, S)=R for S > B, (18) 

where B is the knock-out barrier level and R is the rebate paid in the case of knock-out. In 
terms of the time stepping algorithm: 

= FFT _1 [FFT[u n ] • e* At ] • H B , (19) 



v 



where Hb(x) = l{ a; <in(B/s , (o))} + R " l{z>ln(B/5(o))} ■ By noting that the Fourier transform of a 
product of two functions is the convolution of their respective Fourier transforms, it is possible 
to lift the algorithm to the frequency domain entirely. Thus, the time-stepping algorithm can 
be modified to 

ff 1 - 1 = {FFT[H B }}®{v n -e* At }. (20) 



where the convolution of two vectors, denoted by <S>, can be executed efficiently using the FFT 
algorithm. This approach requires only a single FFT operation (since FFT[Hb] is precomputed) 
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per time-step as compared to two operations in the original approach. This is schematically 
represented in the right panel of Figure [TJ 



Numerical experiments show that direct application of (19) results in erratic convergence, 



especially near the barrier. To stabilize the algorithm we enforce constraint (18) numerically 



via the method of images (see e.g. Buchen (1996) ) by truncating the values of V(t, S) at S = B 



and extending it to an odd function, i.e. setting V(t, B + x) = 2R — V(t, B — x) for x > and 
V(t, B) = R. This procedure is repeated at each time step and, without introducing any bias 
into the solution of the equation for S < B, it stabilizes the convergence of the FST algorithm 



(see Table 15 ). 



3.3.4. Shout Options 

Shout options constitute a significantly more complicated but interesting class of problems 
for numerical valuation. Much like American options, the holder of a shout option may exercise 
their right at any time during the contract lifetime. Specifically, the contract allows the holder 
to "shout" and set a new strike level. The shout regions are obtained by solving a dynamic 
programming problem at each time step similar to computing the optimal exercise boundary of 
American options. An n-shout option allows the investor to "shout" and reset the strike price 
n times, requiring a more subtle analysis in pricing of such options. For n-shout options, the 
FST algorithm must be extended to track several prices at once (each price corresponding to the 
number of shouts remaining) . Details of the application of the FST algorithm to these options 



is forthcoming in Jackson, Jaimungal, and Surkov (2007) 



3.4. Numerical Results 

Below we present our pricing results and compare them to the prices found in the literature. 
The option and stock price models are specified below each table. The closed-form price refers 
to the option price calculated by an analytic formula, if available, while quoted price refers to 
the price calculated by the authors of the paper. Note that the quoted result is the most precise 
value given in the source and not the value that the algorithm converges to. If a closed- form 
formula is not available, a very good approximation for European options can be found by 



evaluating the integral form in Carr and Madan (1999) using an adaptive quadrature method. 



This price is referred to as the integral price. 

It is also of great interest to establish the convergence properties of the FST algorithm. Here 



we follow the estimation approach taken by d'Halluin, Forsyth, and Vetzal (2005) and extend it 



to estimate the order as a function of At and Ax independently. First, assume that 

v approx (At, Ax) = v exac t + c t (At) Pt + c x {Ax) Px , 

where q, pt, c x , and p x are constants. Since the algorithm does not require time-stepping to 
value European options, the equation above can be simplified to depend only on Ax, 

Vapprox(Ax) — V exac t -\- C x ( Ax^ x . 
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N 



Value 



Change log 2 Ratio CPU-Time 



2048 

4096 

8192 

16384 

32768 



18.00329705 
18.00354600 
18.00360820 
18.00362375 
18.00362764 



0.000249 
0.000062 
0.000016 
0.000004 



2.0008 
2.0004 
2.0002 



0.002712 
0.003550 
0.007105 
0.014658 
0.026558 



Table 2 

Option: 
0.15, A = 

Price: 18.0034 Source: Andersen and Andreasen (2000) 



European put S - 
0.1, ft = -1.08, a 



100.0, K = 100.0, 
0.4, r = 0.05, q = 



T = 10; Model: Merton jump-diffusion a = 
0.02; Closed-Form Price: 18.003629 Quoted 



We use estimates of the option price v approx on successively finer grids in space to establish the 
rate of convergence via 

1 | ^approx (Ax) Vapprox (Ax/2) | , , 

Px=log 2 ] , . _ TT .... . (21) 

\Vapprox\'-^X / Z j Vapprox\t-±% I ^)\ 

Here, the absolute changes in the numerator and the denominator are given in the table un- 
der the column "Change" , while the estimated rate of convergence is given under the column 
"log 2 Ratio". 

For European options under various processes (see Tableland Tables [7] - [9] in Appendix |A|) 
we find that the FST algorithm is order 2 in the space variable. For path dependent options it 
is also necessary to establish convergence properties of the algorithm in the time variable. By 
holding Ax constant, the error becomes dependent only on At. We assume 

Vapprox (At) — V eX act Ct(At) 

with pt estimated by 

p — \ „ I 'approx (At) - v approx (At/2)\ (22) 

| Vapprox 

(At/2) - 

Vapprox (At/4)|- v ' 

Results of estimating pt are presented in Appendix [B] and overwhelmingly suggest that the FST 
algorithm is order 1 in the time dimension. Estimating pt and p x independently leads us to 
believe that a good estimate for the error in the FST algorithm is 

Vapprox (At, Ax) — V eX act ^(A)^, 

where At = 0(A 2 ) and Ax = O(A). Again p can be estimated by computing the log ratio of 
changes in v ap prox as At is reduced by the square of the relative reduction in Ax, 

^ I Vapprox (At, Ax) v C[ pp r . 0:r (At/4, Ax/2) I (23) 

P ~ ° g2 \v approx (At/4, Ax/2) - Vapprox (At/16, Ax/4) I ' 
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N 


M 


Value 


Change 


log 2 Ratio 


CPU-Time 


2048 


128 


9.22186786 






0.031823 


4096 


512 


9.22447378 


0.002606 




0.169965 


8192 


2048 


9.22520088 


0.000727 


1.8416 


1.756965 


16384 


8192 


9.22538212 


0.000181 


2.0042 


13.571508 


32768 


32768 


9.22542570 


0.000044 


2.0564 


156.621382 



Table 3 

Option: American put S = 90.0, K - 
191.2, Y = 1.0102, r = 0.1; Quoted Price: 9.2185 Source: |Forsyth, Wan, and Wang (2006) 



1.0, T = 0.25; Model: CGMY C = 0.42, G = 4.37, M 



N 


M 


Value 


Change 


log 2 Ratio 


CPU-Time 


2048 


128 


0.25348926 






0.016819 


4096 


512 


0.25401054 


0.000521 




0.181466 


8192 


2048 


0.25414993 


0.000139 


1.9030 


1.297593 


16384 


8192 


0.25418356 


0.000034 


2.0514 


10.794937 


32768 


32768 


0.25419310 


0.000010 


1.8178 


123.635629 



Table 4 

Option: Up-and-Out Barrier Call S = 100.0, if = 100.0,5 = 110, T = 1.0; M odel: Black- 
Scholes-Merton a = 0.15, r = 0.05, q = 0.02; Closed-Form Price: 0.2541963 Source: |Hull (2005) 



N M Value Change log 2 Ratio CPU-Time 

2048 128 35.73104311 0.719093 

4096 512 35.88365480 0.152612 4.285660 

8192 2048 35.92030459 0.036650 2.0580 31.988542 

16384 8192 35.93057212 0.010268 1.8357 286.700219 

Table 5 

Option: Shout S = 100.0, T = 10.0, n = 5; Model: CGMY C = 1.0, G = 5.0, M = 5.0, Y = 
0.5, r = 0.1. 



The results presented in Tables [3] -|5] and in Tables 10 - 12 in Appendix |A| support the conclusion 
that the error in the FST algorithm behaves as follows 



J approx 



(At, Ax) = v exact + c t At + c x (Ax) 



(24) 



Our implementation of the FST algorithm uses the FFTW library developed by Frigo and 



Johnson (2005) to compute the FFTs. The numerical experiments were written in C++ and 
performed on a Pentium 4 1.50 GHz machine. 



4. Extension to Multi- Asset Problems 



Pricing of contingent claims that depend on two or more assets is a complicated task. Typical 
approaches include Monte- Carlo simulations, solving multi-dimensional PDEs and FFT based 
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approaches. Although very flexible, Monte-Carlo methods can suffer from slow convergence 
when jumps are incorporated into the model. Solving a multi-dimensional PDE is more efficient; 
however, when jumps are introduced, the resulting PIDEs are difficult to handle due to the non- 
local integral term. Based on the FFT approach of Carr and Madan (1999) |Dempster and 



Hong (2000) develop a method for pricing spread options - options on two underlying stocks. 
Unfortunately, their method is only applicable to European options and does not readily extend 
to path-dependent options. 

The computational efficiency of the FFT is not restricted to a single dimension; as such, 
the multi-dimensional FFT can be leveraged to efficiently compute multi-dimensional DFT. 
Therefore, it seems natural to generalize our method to a multi-dimensional FST algorithm. 
The advantages of our approach over previous approaches is its ease in handling any Levy 
process and its ease of application to path-dependent options. In this section, we present the 
main ideas of the derivation and conclude with pricing results for two- asset European, American 
and barrier options. The mathematical basis of our method follows through from Section [2] 

Let V(t,S(t)) denote the price at time t of an option, written on a vector of d underlying 
price indices S(t), whose components are Sj(t), with a T-maturity payoff of (p(S(T)). The value 
of the option is the discounted expectation under the risk-neutral measure Q 

V(t, S(t)) = Ef [e- p C r -*> <p(S(T)j\ . (25) 

The discount-adjusted and log-transformed price process v(t, X(t)) := e r ( T ~^V(t, S(0)e x ^) 
is a Q-martingale. We assume that the underlying indices S(i) follow a ci-dimensional expo- 
nential Levy process with a characteristic triplet (~f,C,v), where 7 represents the vector of 
unadjusted-drifts, C represents the variance-covariance matrix of the diffusions, and v is the 
multi-dimensional Levy density. As before, X admits a canonical Levy-Ito decomposition 

X(t) = 7 t + W(t) + 3 l (t) + lim J e (t) , (26) 

e\0 

At) =11 YKdyxds), (27) 

J0 J\y\>l 

F(t) = f f y Hdy x ds) - u(dy x ds)] . (28) 

JO Je<\y\<l 

Applying the zero-drift condition on v together with its boundary condition at maturity leads 
to the PIDE 

j(a + £)v=o, 

\ k(T,x) = v (S(0)e*), 
where C is the infinitesimal generator of the multi-dimensional Levy process and acts on twice 
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differential) le functions /(x) as follows 



£/(x) = (7 • d x + \d x ■ C ■ d x ) /(x)+ / (/(x + y) - /(x) - y • d x /(x) l, y|<1 ) !/(dy) (30) 

JR™/{0} 



A d-parameter family of ODEs is obtained by applying Fourier transforms to both sides of (29) 
8^[«](i,w) + $(«)5[«](t,w) = 0. (31) 

Here, 

*M =7'U + lwC«+ / (e iu " y - 1 - iy • w l, y|<1 ) (32) 

is the characteristic exponent of the d-dimensional Levy process. The system of ODEs admits 
the simple analytic solution 

T[v](t 1 ,u)=T[v](t 2 ,w) e^-* 1 ^ (33) 

resulting in a natural generalization of the FST time-stepping method to arbitrary dimensions 

FFT~ 1 [FFT[u n ] • e* A *] . (34) 



n-l 



Of course, this algorithm could be derived without going through the PIDE, and instead relying 
solely on the martingale property of v. 

If jumps in individual assets are uncorrelated, then the multi-dimensional Levy density factors 
v{dy x ds) = V\(dy\) . . . v n {dy n ) ds , leading to a factorization of the characteristic function. 
However, allowing for correlation between the jumps in the assets is just as straightforward as 



long as the integral appearing in (32 ) can be computed analytically in closed form. In an extreme 
case, there may be a single market wide jump risk factor together with idiosyncratic jump risks. 
In the next subsection we illustrate two applications of the FST algorithm in two-dimensions. 



4.1. Spread Options 

An interesting class of multi-asset options are spread options - the option to exchange /3-units 
of one asset for a-units of another asset. These options can be viewed as options on the difference 
(or spread) of two stock prices with terminal payoff 

p(Si(T), S 2 (T)) = max(aS 2 (T) - 0S\(T) - K, 0) . (35) 

Spread options do not admit an analytic closed-form solution even for the Black-Scholes-Merton 



model. For a detailed discussion of spread options and various approximations see Carmona and 



Durrleman (2003) |Dempster and Hong (2000) present an FFT-based approach to valuation of 



spread options. Their approach involves breaking the region in which the option is in-the-money 
into a series of rectangular approximations. Unfortunately, they only apply their method to a 
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Figure 2. Price surface of a European spread call with parameters T = 1, a = 1.0, S\ = 96, (3 
1.0, 6*2 = 100, K = 0. The stock price process is Merton jump-diffusion r = 0.1, o\ = 0.1, gi 
0.05, Ai = 0.25, fix = -0.13, d\ = 0.37, a 2 = 0.2, q 2 = 0.05, A 2 = 0.5, fi 2 = 0.29, <f 2 = 0.41, p 
0.5 



pure diffusion model with stochastic volatility and it seems difficult to extend this method to 
the Bermudan or barrier cases. 

For our numerical experiments, we assume a joint jump-diffusion with Merton-like jumps and 



compare our results with the Kirk (1995) approximation and its extension for jump-diffusions 



found in Carmona and Durrleman (2003) In this case, the Levy density factors with Vi{dy) = 
(Ai/ \j2irdf) exp{—(y— pi) 2 /2af}dy and the diffusive volatilities are 0{ with correlation p. Hence, 



i(px 



a 



2 

Mwi + i(p 2 - ^f) 



0J 2 



+ Ai ( e 'Ai^i-5-^i/2 



2' 2 
1) + A 2 (e*' 



paia 2 iOiio 2 
1) 



(36) 



where the drift is fixed by risk-neutrality to be pi = r — Aj (e^ +CT i 2 / 2 — 1). In Appendix [c] 
we present numerical results for European and American options. Of course, the method is 
applicable to barrier spread options as well. 
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4.2. Catastrophe Options 

Catastrophe options have become more important in recent years, yet there is very little 
published on the subject. These options pay the holder a function of total losses and the 
company's equity value. As a result, it is important to jointly model losses and equity, especially 



since large losses will cause significant drops in share value. Cox, Fairchild, and Pedersen (2004) 



introduce a simple model which counts the number of losses paying no attention to loss size. 



Jaimungal and Wang (2006) extend the model to incorporate random losses as well as stochastic 



interest rates. However, both works assume that the option is European, while the existing 
contracts are in fact of American type. Here, we illustrate how the FST can be used to price 
the early exercise premium. In the constant interest case, |Jaimungal and Wang (2006) make the 
following assumption on share value S(t) and losses L(t) 



S(t) = 5(0)exp{J(t) + 1 t + aW t } , (37) 

J(t) = -aL(t), (38) 

N(t) 

L(t) = ( 39 ) 



n=l 



where N(t) is a compound Poisson process with activity rate A, li are i.i.d. random variables 
with probability density //,(/) having support on M + . Notice that it is the presence of losses 
which drives the jumps in the price process and not an independent jump process. 

In this case, the 2-dimensional Levy density is v{dy\ x dy<i) = fL(y2)$(yi + ciyi) dy\dy2 
resulting in the characteristic function 



^(wi.wa) = i 7 wi - \o 2 u{ + [ (e i( -- auJ1+ ^ y - 1 - iy (-a^ + cu 2 ) t lyl<1 ) dy . (40) 

J —oo 

The risk-neutral drift is 7 = r — Vl/(— i,0). With this characteristic function a slight modifica- 
tion of the American styled options algorithm leads to an efficient pricing mechanism. Namely, 
the exercise policy is now chosen at each point in the (S(t), L{t)) plane independently. If L(t) was 
not a separate observable, as it is in the usual jump-diffusion model case, then the exercise policy 
would be independent of L{t). We leave the application of the FST method to catastrophe op- 
tions such as catastrophe equity put options (with payoff cp(S(T), Lt) = tL T >L t +u(K — S(T)) + ) 
for future work. 



5. Regime Switching Models 



Regime switching models can be traced back to the early work of Lindgren (1978) and ever 
since the seminal work of Hamilton (1989 1990) have become a very popular approach to 
incorporate non-stationary behaviour into an otherwise stationary model. The essential idea 
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is to assume that the world switches between states representing, for example, moderate, low 
and high volatilities regimes. Although popular for describing time-series, little work has been 
carried out in terms of option valuation. Two-state European options in log-normal models 



were studied in Naik (1993); while European options in a two-state VG model were studied 



by Konikov and Madan (2000) |Albanese, Jaimungal, and Rubisov (2003) derive closed form 



results for barrier and European options, and semi-closed form formulae for American options, 
in a special class of two-state VG models. Here, we demonstrate that the FST algorithm can 
easily incorporate path-dependent options, such as barrier and American options, with multiple 
regimes and multiple assets in computational time proportional to the number of regimes. 

Let K := {1, . . . , K} denote the possible hidden states of the world, and let Z{t) G K denote 
the prevailing state of the world at time t. We will assume that Z(t) is driven by a continuous 
time Markov chain with generator A, i.e. the transition probability from state k at time t\ to 
state I at time t 2 is P^ 2 := Q(Z(t 2 ) = l\Z(tx) = k) = (exp{(t 2 - h)A}) kl . The real matrix A 
satisfies the usual requirements: An = — Ylk^i ^ik and A k i > VTc ^ I. Given that Z(t) = k, we 
assume that the joint stock price process S(i) follows a (i-dimensional exponential Levy process 
with Levy triple (■y^ k \ C^ k \ v^). The drift vectors of each state are assumed prefixed at their 
risk-neutral levels of jj = r — ty( k \—ilj), where \l/( fc )(u>) denotes the characteristic exponent of 
the respective Levy processes and 1,- is the vector with zeroes everywhere except a single entry 
of 1 at dimension j. This modeling assumption can succinctly be written cZX(i) = dXS z ^'\t), 
where ~X.( k \t) is the k-th d-dimensional Levy process and the price processes are obtained by 



exponentiation component wise: Sj(t) = Sj(0)exp{Xj(t)}. Chourdakis (2005) investigates the 



d = 1 version of this framework and derives the characteristic function of the terminal stock 
price. 

The author calculates European option prices via FFT methods; however, then resorts to 
numerical integration for the valuation of path-dependent options. We take a slightly different 
approach, and make use of a generalization of the FST algorithm which allows path-dependent 
options based on the regime switching models to be valued efficiently. 

Under the above assumptions, let u(X(t), Z(t), t) denote the discounted-adjusted and log- 
transformed price at time t conditional on the state Z(t) and spot levels X(t). It is not difficult 
to show that European option prices satisfy the following system of PIDEs: 

{ d t +(A kk + C^)v(x,k,t) + , £ jjik A jk v(x,j ) t) =0, 

\ v(x,k,T) = ^(S(0)e x ) , 

for every k E K. Here, represents the infinitesimal generator of the k-th d-dimensional 
Levy process. It is possible in principle to apply any of the usual finite-difference schemes to 
this system of PIDEs to solve the problem. However, as discussed earlier, this is quite difficult 
due to the non-local integral terms and especially so for multi-dimensional problems. Instead, 
we develop an FST algorithm. First discretize the continuous time Markov chain in the usual 
manner by partitioning time into steps of size At and assuming Z(t) is held constant on time 
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Figure 3. This diagram depicts the effect of introduction of regime switching on option price 
of an American option. Option: American put S = 100.0, -fT = 100.0, T = 1.0; Model: Merton 
jump-diffusion a = 0.15, fl = —0.5,5" = 0.45, r = 0.05 Regime Switching: parameter states 
A € [0.3,0.5,0.7], initial state probabilities q = [0.2,0.3,0.5], Markov chain generator A = 
[-0.4, 0.3, 0.1; 0.1, -0.5, 0.4; 0.05, 0.15, -0.2] 



intervals (i n _ i,t n ] (with t n = nAt), for n £ N, with transition probabilities 
Pkl : = 



A kl At , k^l , 
1 + An At , otherwise . 



(42) 



Then, by the martingale property of v and the law of iterated expectations we have 
v(X(t n ), Z(t n ),t n ) = E? [v(X.(t n +i), Z(t n+ i),t n+ i)] 



E Q [v(X(t n+1 ),Z(t n+1 ),t n+1 )\F tn V Z(t n+1 )} 



Within the inner expectation, the process X(t) follows a given d-dimensional Levy model, 
resulting in an expectation of the single regime form; consequently, the inner expectation can 
be written 



F[v](u),Z(t n+ i),t n+1 ) exp 



(43) 
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Figure 4. This diagram depicts the effect of introduction of regime switching on exercise bound- 
ary of an American option. Option: American put S = 100.0, K = 100.0, T = 1.0; Model: 
Merton jump-diffusion a = 0.15, fi = —0.5,5" = 0.45, r = 0.05 Regime Switching: parameter 
states A £ [0.3,0.5,0.7], initial state probabilities q = [0.2,0.3,0.5], Markov chain generator 
A = [-0.4, 0.3, 0.1; 0.1, -0.5, 0.4; 0.05, 0.15, -0.2] 



The above representation can be interpreted as the (discount-adjusted) price of a contingent 
claim with payoff u(X(i n +i), Z(t n +i), i n +l) assuming that the multi-dimensional price process is 
following c?-dimensional Levy model Z(t n +i). Finally, due to the linearity of the inverse Fourier 
transform, the outer expectation can be computed to obtain the simple iterative scheme 

K 

v(X(t n ),j,t n ) = ^Pjk ?~ 

k=l 

At each time step, the algorithm therefore requires storing K prices. These K prices are then 
integrated backwards in time by the FST algorithm, then weighted according to the transition 
probabilities. If there are exercise decisions to be made, these must be made after averaging. 
In this manner, the price of the option today in all K states will be known. Since a trader will 
likely not know for certain what state the world is currently in, the price will not simply be 
one of these prices; rather, the trader must decide, exogenously, on a probability qj~ that the 



F[v](u,k,t n+1 ) exp{Ai*( fc V)} (X(i n )) . 



(44) 
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N 


M 


Value 


Change 


log 2 Ratio 


CPU- Time 


2048 


128 


33.76666279 






0.104524 


4096 


512 


33.77807492 


0.011412 




0.693101 


8192 


2048 


33.78118573 


0.003111 


1.8752 


6.960128 


16384 


8192 


33.78168754 


0.000502 


2.6321 


64.977226 


32768 


32768 


33.78181086 


0.000123 


2.0248 


656.755101 



Table 6 

Option: American put S = 1369.41, if = 1200.0, T = 0.56164; Model: Variance- Gamma 
fj, = —0.22898,cr = 0.20722, r = 0.0541, q = 0.012; Regime Switching: parameter states k £ 
[0.22460,0.49083,0.50215], initial state probabilities q = [0.2,0.3,0.5], Markov chain generator 
A = [-0.4, 0.3, 0.1; 0.1, -0.5, 0.4; 0.05, 0.15, -0.2] 



world is in state k. Once these probabilities are determined, the trader's price for the option is 

££=i?fc«(x(o),fc,o). 

In Table [6] we present the pricing results of an American put option with regime switching. 
The price of the option with regime switching (33.781) is lower than the price of same American 
option with k = 0.50215 and without regime switching (35.524). This is expected since the 
former option switches between periods of low, medium and high subordinator volatility while 
the latter option is always in the state of high subordinator volatility. 

6. Conclusions and Future Work 

We introduced a new method for pricing European and path-dependent options when the 
underlying process(es) follows a regime switching Levy process(es). The method treats the 
integral term and diffusion terms in the pricing PIDE symmetrically, is efficient, and accurate. 

The numerical results attained by the FST method are quite impressive. We succeeded in re- 
producing the pricing results obtained by various authors. For European options, our algorithm 
is preferable to PIDE finite-difference based methods since it does not require a time-stepping 
procedure. Our algorithm is more appealing than the usual Fourier transform methods since 
the analytic expression for the Fourier transform of the option payoff is not required. Thus, 
the FST method is the method of choice for European options, especially with non-standard 
payoffs. We showed that Barrier options are priced accurately by comparing with analytical 
results in the BSM model. Furthermore, we compare our results for Barrier options under jump 
models with those of other authors and find excellent agreement. Additionally, we matched 
results obtained by other published methods for American options. This confirms that the FST 
is a robust method for pricing options with excellent precision properties. Moreover, the FST 
method provides a generic framework for option pricing under any stock price process, such 
as Brownian motion, jump-diffusion or exponential-Levy. We also demonstrated how the FST 
method easily extends to the multi-dimensional case and can incorporate regime switching. 

At this stage of development, the FST algorithm is first order in time. It should be noted, 
however, that second order methods in time typically involve a multi-step or iterative refinement 
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of the solution at each time step. It is of great interest to investigate various approaches to im- 
prove the FST time order. For American options, this may involve an iterative approach. Barrier 
options are more challenging than American options due to the discontinuity in the option value 
created by the knock-out / knock-in provisions. In finite-difference methods, this problem is 
addressed by refining the mesh around the point of discontinuity. FFT algorithms, however, are 
designed for uniformly spaced points, making mesh refinement difficult to implement. Fortu- 



nately, FFTs using unequal spaced data (NFFT) have been developed recently. Potts, Steidl. 



and Tasche (2000) provide an overview of existing NDFT methods and develop a method of 



their own. An interesting area of research would be the use of NFFT algorithms in the FST 
method to improve its time-stepping order. Other areas of potential research include: compu- 
tation of American exercise boundary, efficient computation of the Greeks, processes parameter 
calibration, and stochastic optimal control problems. 
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A. Pricing Results 



N 



Value 



Change log 2 Ratio CPU- Time 



2048 

4096 

8192 

16384 

32768 



0.04261423 
0.04263998 
0.04264641 
0.04264801 
0.04264841 



0.000026 
0.000006 
0.000002 
0.000000 



2.0018 
2.0010 
2.0011 



0.002460 
0.005240 
0.009759 
0.019089 
0.038259 



Table 7 

Option: European call S = 1.0, K = 1.0, T = 0.2; Model: Kou jump-diffusion u = 0.2, A = 
0.2, p = 0.5, 77_ = 3,t?+ = 2,r = 0.0; Integral Price: 0.0426478 Quoted Price: 0.0426761 Source: 
Almendral and Oosterlee (2005) | 



N 



Value 



Change log 2 Ratio CPU- Time 



2048 

4096 

8192 

16384 

32768 



7.49983308 
7.50061473 
7.50091004 
7.50098386 
7.50100232 



0.000782 
0.000295 
0.000074 
0.000018 



1.4043 
2.0000 
2.0002 



0.002743 
0.005478 
0.011147 
0.024923 
0.043864 



Table 8 

Option: European call S = 100.0, K = 100.0, T 
-0.28113,ct = 0.19071, k = 0.49083, r = 0.0549, q = 



= 0.46575; Model: Variance-Gamma fj, = 
0.011; Integral Price: 7.50100847 ^owrce: 



Hirsa and Madan (2004) 



N 



Value 



Change log 2 Ratio CPU-Time 



2048 

4096 

8192 

16384 

32768 



0.10295883 
0.10296489 
0.10296640 
0.10296678 
0.10296687 



0.000006 
0.000002 
0.000000 
0.000000 



2.0009 
2.0004 
2.0002 



0.005580 
0.009415 
0.018672 
0.031994 
0.056374 



Table 9 

Option: European put S 
5.0, Y = 0.5, r = 0.1; Integral Price: 0.10296691 Quoted Price: 0.1029669 Source: |Almendral 



1.0, K = 1.0, T = 1.0; Model: CGMY C = 1.0, G = 5.0, M 



and Oosterlee (2007) 
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N 


M 


Value 


Change 


log 2 Ratio 


CPU-Time 


2048 


128 


3.23945333 






0.020389 


4096 


512 


3.24080513 


0.001352 




0.169292 


8192 


2048 


3.24114185 


0.000337 


2.0053 


1.736765 


16384 


8192 


3.24122597 


0.000084 


2.0011 


14.570725 


32768 


32768 


3.24124692 


0.000021 


2.0050 


150.485790 



Table 10 

Option: American put S = 100.0, K = 100.0, T = 0.25; Model: Merton jump- diffusion a = 
0.15, A = 0.1, fx = -0.9, a = 0.45, r = 0.05; Quoted Price: 3.2412435 Source: |d'Halluin, Forsyth, 



and Labahn (2003) 



N 


M 


Value 


Change 


log 2 Ratio 


CPU-Time 


2048 


128 


35.50814417 






0.021277 


4096 


512 


35.51990755 


0.011763 




0.162960 


8192 


2048 


35.52355311 


0.003646 


1.6901 


1.646470 


16384 


8192 


35.52432105 


0.000768 


2.2471 


15.200406 


32768 


32768 


35.52449938 


0.000178 


2.1065 


142.993274 



Table 11 

Option: American put S = 1369.41, K = 1200.0, T = 0.56164; Model: Variance-Gamma 
H = -0.22898, a = 0.207 22, k = 0.50215, r = 0.0541, q = 0.012; Quoted Price: 35.5301 Source: 
Hirsa and Madan (2004) 



N 


M 


Value 


Change 


log 2 Ratio 


CPU-Time 


2048 


128 


8.89073604 






0.026034 


4096 


512 


8.89040553 


0.000331 




0.188495 


8192 


2048 


8.89031967 


0.000086 


1.9446 


1.261334 


16384 


8192 


8.89029787 


0.000022 


1.9781 


10.623886 


32768 


32768 


8.89029207 


0.000006 


1.9105 


126.700997 



Table 12 

Option: Down-and-Out Barrier Call S = 100.0, K = 110.0,5 = 85, T = 1.0; Model: Merton 
jump-diffusion a = 0.25, A = 2.0, jl = 0.0,5- = 0.1, r = 0.05; Quoted Price: 9.013 Source: 
Metwally and Atiya (2003) 
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B. Time-Stepping Results 



N 



Value 



Change log 2 Ratio CPU-Time 



4096 
4096 
4096 
4096 



512 
1024 
2048 
4096 



35.51990755 
35.52162452 
35.52246665 
35.52288380 



0.001717 
0.000842 
0.000417 



1.0277 
1.0135 



8192 
8192 
8192 
8192 



512 
1024 
2048 
4096 



35.52073332 
35.52261232 
35.52355311 
35.52402453 



0.001879 
0.000941 
0.000471 



0.9980 
0.9969 



16384 
16384 
16384 
16384 



512 
1024 
2048 
4096 



35.52079913 
35.52267988 
35.52361898 
35.52408731 



0.001881 
0.000939 
0.000468 



1.0020 
1.0038 



Table 13 

Option: American put S = 1369.41, if = 1200.0, T 
H = -0.22898,ct = 0.20722, k = 0.50215, r = 0.0541, q = 



0.147402 
0.495313 
0.651053 
1.427344 



0.403792 
0.852099 
1.612530 
3.503661 



0.989980 
1.735797 
3.995848 
7.307692 



= 0.56164; Model: Variance- Gamma 
0.012; Quoted Price: 35.5301 Source: 



Hirsa and Madan (2004) 



N 

4096 
4096 
4096 
4096 



M 

512 
1024 
2048 
4096 



Value Change log 2 Ratio CPU-Time 



9.22447378 
9.22493756 
9.22517293 
9.22529067 



0.000464 
0.000235 
0.000118 



0.9785 
0.9993 



8192 
8192 
8192 
8192 



512 
1024 
2048 
4096 



9.22451151 
9.22497177 
9.22520088 
9.22531447 



0.000460 
0.000229 
0.000114 



1.0064 
1.0122 



16384 
16384 
16384 
16384 



512 
1024 
2048 
4096 



9.22451347 
9.22497714 
9.22520900 
9.22532453 



0.000464 
0.000232 
0.000116 



0.9999 
1.0050 



0.347249 
0.874211 
1.002758 
1.798958 



0.689149 
1.288446 
1.843247 
4.063315 



1.172726 
2.521093 
4.401092 
9.088098 



Table 14 

Option: American put S = 90.0, K = 98.0, T = 0.25; Model: CGMY C = 0.42, G = 4.37, M 



191.2, Y = 1.0102, r = 0.1; Quoted Price: 9.2185 Source: Forsyth, Wan, and Wang (2006) 
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N 

4096 
4096 
4096 
4096 



8192 
8192 
8192 
8192 



16384 
16384 
16384 
16384 



M 

512 
1024 
2048 
4096 



Value Change log 2 Ratio CPU-Time 



0.25401054 
0.25409348 
0.25413559 
0.25415687 



0.000083 
0.000042 
0.000021 



0.9779 
0.9844 



512 
1024 
2048 
4096 



0.25402488 
0.25410782 
0.25414993 
0.25417121 



0.000083 
0.000042 
0.000021 



0.9779 
0.9844 



512 
1024 
2048 
4096 



0.25402651 
0.25410944 
0.25415155 
0.25417284 



0.000083 
0.000042 
0.000021 



0.9779 
0.9844 



0.168113 
0.310622 
0.717033 
1.426751 



0.466672 
0.985679 
1.661753 
3.377130 



1.070065 
2.567016 
4.697512 
7.583904 



Table 15 

Option: Up-and-Out Barrier Call S = 100.0, K = 100.0,5 = 110.0, T = 1.0; Model: Black- 
Scholes-Merton a = 0.15, r = 0.05, q = 0.02; Closed-Form Price: 0.2541963 Source: |Hull (2005) 



N M Value Change log 2 Ratio CPU-Time 

4096 512 8.89040553 0.170580 

4096 1024 8.89035362 0.000052 0.412546 

4096 2048 8.89032768 0.000026 1.0006 0.890808 

4096 4096 8.89031471 0.000013 1.0003 1.822582 

8192 512 8.89039752 0.432691 

8192 1024 8.89034561 0.000052 1.296842 

8192 2048 8.89031967 0.000026 1.0006 1.600545 

8192 4096 8.89030670 0.000013 1.0003 3.302225 

16384 512 8.89039517 1.257726 

16384 1024 8.89034327 0.000052 2.335755 

16384 2048 8.89031732 0.000026 1.0005 4.781386 

16384 4096 8.89030435 0.000013 1.0004 9.956966 

Table 16 

Option: Down-and-Out Barrier Call S = 100.0, K = 110.0,5 = 85.0, T = 1.0; Model: Merton 

jump-diffusion a = 0.25, A = 2.0, fl = 0.0,5- = 0.1, r = 0.05; Quoted Price: 9.013 Source: 
Metwally and Atiya (2003) 
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C. Multi- Asset Pricing Results 



N Value Change log 2 Ratio CPU-Time 



512 

1024 

2048 

4096 

8192 



7.55540525 
7.53923270 
7.54214361 
7.54233725 
7.54232390 



0.016173 
0.002911 
0.000194 
0.000013 



2.4740 
3.9100 
3.8581 



0.450042 
1.944840 
8.654331 
37.481155 
159.032330 



Table 17 

Option: Spread call a = 1.0, Si = 96.0,0 = 1.0, S 2 = 100.0, K = 
Scholes-Merton a x = 0.1, gi = 0.05, o 2 = 0.2, q 2 = 0.05, p = 0.5, r 



2.0, T = 1.0; Model: Black- 
= 0.1; Kirk's Formula Price: 



7.54232193 Quoted Price: 7.542242 Source: Dempster and Hong (2000) 



N 



Value 



Change log 2 Ratio CPU-Time 



512 

1024 

2048 

4096 

8192 



15.03639950 
15.02776432 
15.02919574 
15.02924971 
15.02924214 



0.008635 
0.001431 
0.000054 
0.000008 



2.5928 
4.7293 
2.8345 



0.880245 
2.821585 
11.919293 
48.371978 
209.806361 



Table 18 

Option: Spread call Si = 96.0, S 2 = 100.0, K = 
<T x = 0.1,gi = 0.05, Ai = 0.25, fix = -0.13, a x 



2.0, T = 1.0; Model: Merton jump-diffusion 
= 0.37, 0-2 = 0.2,^2 = 0.05, A 2 = 0.5, fi 2 = 



0.11,(72 = 0.41, p = 0.5, r = 0.1; Kirk's Formula Price: 15.03001533 



N M Value Change log 2 Ratio CPU- Time 

512 64 5.61098208 3.009590 
1024 256 5.61273425 0.001752 51.875703 
2048 1024 5.61216093 0.000573 1.6117 1373.959488 
4096 4096 5.61199872 0.000162 1.8215 30498.836319 

Table 19 1 

Option: American spread put a = 1.0, Si = 96.0, = 1.0, S 2 = 100.0, K = 2.0, T = 1.0; Model: 
Black-Scholes-Merton a x = 0.1, q x = 0.05, a 2 = 0.2, q 2 = 0.05, p = 0.5, r = 0.1; 
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